Introducción
En esta práctica vamos a tratar con datos de defunciones por enfermedad isquémica en Aragón, separados por municipio. Nuestro objetivo va a ser modelizar el riesgo de muerte a causa de esta enfermedad, y realizar una predicción del riesgo “real” de muerte para cada municipio.
1. Lectura y preparación de los datos
# Librerías
# install.packages("rmdformats") # plantillas para reporte HTML
library(dplyr) # operaciones con tablas
library(sf) # lectura y tratamiento de datos geoespaciales
library(spdep) # lectura y tratamiento de datos geoespaciales
library(R2WinBUGS) # modelos bayesianos, llama a WinBUGS
library(INLA) # modelos bayesianos utilizando INLA
library(ggplot2) # representación gráfica
theme_set(theme_void()) # cambiar estética de todos los plots de ggplot
library(gridExtra) # extensión de ggplot para plotear varias gráficas a la vez
library(leaflet) # representación de mapas interactivos
Comenzamos cargando las librerías necesarias para nuestra práctica y leyendo los archivos de datos:
# Cargamos los datos con info por municipio
load(file.path("datos", "Aragon.Rdata"))
datos <- Aragon.df
# Leemos los datos georeferenciamos y los unimos con los datos cargados
mapa <- read_sf(file.path("datos", "aragon.shp")) %>%
left_join(datos, by = "CODMUNI") %>%
mutate(RME = O / E)
Para cada municipio tenemos el número de muertes observadas por enfermedad isquémica en hombres, así como el número esperado de muertes, calculado como el número de muertes que le correspondería a cada municipio dada su población, multiplicando la población del municipio por el número promedio de defunciones por habitante de Aragón. Como una primera aproximación, vamos a representar los observados, los esperados y el riesgo de mortalidad estandarizado en un mapa:
Se puede apreciar en los mapas del número de casos observados y esperados que el número de fallecimientos es muy dependiente de la población de cada municipio, y que si no lo ajustamos de alguna forma lo único que vamos a observar en el mapa es un reflejo de la población de cada municipio. Al dividir observados entre esperados obtenemos la Razón de Mortalidad Estandarizada (RME), que ya es un indicador más preciso del riesgo real de fallecimiento. Un RME mayor a 1 indica un exceso de riesgo, por lo que el mapa en el que representamos el RME es útil para conocer las zonas de mayor riesgo. El problema es que muchos municipios de Aragón poseen muy pocos habitantes, por lo que una desviación muy pequeña entre observados y esperados causa variaciones muy grandes en el RME, teniendo zonas con un riesgo aparente de 0, y zonas con un exceso de riesgo enorme, cuando realmente sólo hay uno o dos fallecimientos de diferencia entre observados y esperados. Es por esto que vamos a ajustar un modelo que nos suavice este riesgo, teniendo el cuenta los casos observados en los municipios vecinos, para tratar de encontrar una estimación del riesgo que se asemeje más a la realidad.
2. Ajuste de modelo de Besag, York y Mollié
Vamos a implementar el modelo de Besag, York y Mollié (BYS), tanto en WinBUGS como en INLA. Este es un modelo jerárquico que combina un modelo poisson gamma con un efecto autorregresivo que confiere un efecto espacial, al tener en cuenta las observaciones en las localidades colindantes.
2.1. Modelo en WinBUGS
Ajustamos primero el modelo en WinBUGS, que se basa en el uso de simulaciones mediante métodos MCMC (Marcov chain Monte Carlo). Para ello, nos aseguramos de que el número de iteraciones, el burn-in (número de muestras simuladas iniciales que eliminamos) y el thin (se guarda una muestra cada k iteraciones) son los suficientes para que las cadenas converjan correctamente. El código utilizado para ajustar el modelo es el siguiente:
# Generamos estructura de vecindad
vecinos <- poly2nb(mapa) %>%
nb2WB()
# Modelo BYS en WinBUGS
set.seed(1234)
modeloBugs <- function(){
# verosimilitud
for(i in 1:n){
O[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + m + het[i] + sp[i]
het[i] ~ dnorm(0, prechet)
R[i] <- exp(m + het[i] + sp[i]) # RMS
}
for (i in 1:n) {
pR[i] <- step(R[i]-1) # probabilidad de que riesgo superior a 1
}
# distribuciones a priori
sp[1:n] ~ car.normal(adj[], w[], num[], precsp)
m ~ dflat()
prechet <- pow(sdhet, -2)
precsp <- pow(sdsp, -2)
sdhet ~ dunif(0, 10)
sdsp ~ dunif(0, 10)
# distribución predictiva
for (i in 1:n){
yPred[i] ~ dpois(mu[i])
res[i] <- yPred[i] - O[i]
pRes[i] <- step(res[i])
pY[i] <- step(yPred[i])
}
}
datos<-list(O = mapa$O,
E = mapa$E,
n = nrow(mapa),
adj = vecinos$adj,
w = vecinos$weights,
num = vecinos$num)
iniciales<-function(){
list(m = rnorm(1),
het= rnorm(nrow(mapa),0,0.1) ,
sp = rnorm(nrow(mapa),0,0.1),
sdhet=runif(1,0,1),
sdsp=runif(1,0,1),
yPred = rpois(nrow(mapa),mean(mapa$O)))
}
param <- c("R","m", "het","sp","sdhet","sdsp","pR","pRes","yPred","pY")
resultados <- bugs(
data = datos,
inits = iniciales,
parameters = param,
model = modeloBugs,
n.burnin = 3000,
n.thin = 25,
n.iter = 10000
# debug=TRUE
)
Representamos en un mapa el Riesgo de Mortalidad Suavizado (RMS) para cada municipio, así como la probabilidad estimada de que el riesgo sea superior a 1, que no es sino la proporción de iteraciones en las que el valor estimado por el modelo supera la unidad.
Observamos que el riesgo es inferior en el norte y el sur de la comunidad con respecto al centro, siendo este especialmente alto en el oeste, en la zona de las conocidas “cinco villas” como Tauste, Sádaba y Ejea de los Caballeros. Procedemos a ajustar este mismo modelo utilizando INLA para ver si obtenemos resultados similares. Comparamos las predicciones de observados frente a los valores observados reales:
2.2. Modelo en INLA
# Generamos estructura de vecindad y la guardamos en un archivo
# poly2nb(mapa) %>%
# nb2INLA(file = file.path("datos", "aragon.graph"))
# Modelo BYS en INLA
H <- inla.read.graph(filename = file.path("datos", "aragon.graph"))
S <- U <- seq(1, 729)
mapa <- mapa %>%
mutate(S = S,
U = U)
formula <-
O ~ 1 +
f(S,
model = "besag",
graph = H,
scale.model = TRUE,
hyper = list(
prec = list(prior = "loggamma", param = c(1, 0.001)))
) +
f(U,
model = "iid",
hyper = list(prec = list(prior = "loggamma", param = c(1, 0.001)))
)
modeloINLA <- inla(
formula = formula,
family = "poisson",
data = mapa,
E = E,
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
control.predictor = list(compute = TRUE, cdf = c(log(1)))
)
Representamos de nuevo en un mapa tanto el RMS como la probabilidad estimada de que el riesgo sea superior a 1.
Como se puede apreciar, obtenemos unos riesgos muy similares a los obtenidos con WinBUGS, y, si bien en este caso estos quedan algo más suavizados, en general el riesgo se distribuye de la misma forma en ambas implementaciones. Comparamos de nuevo las predicciones de observados frente a los valores observados reales:
3. Mapa interactivo con librería leaflet
Vamos a representar en un mapa interactivo todos los riesgos y probabilidades estimados, para tener toda la información obtenida a lo largo de la práctica en un único lugar y poder visualizarla y compararla más eficientemente.
Reproducibilidad
La elaboración de este informe ha sido gestionada mediante el uso de git. Puedes descargar todo el código y los datos en este link, o, si eres usuario de git, utilizando el comando git clone https://github.com/Julian-Guillo/Disease-Mapping-Aragon.
Con respecto al código, este ha sido extensamente comentado, los nombres de las variables son representativos de su contenido y siguen la convención camelCase. Para cualquier corrección o duda por favor abrir un pull request o contactar conmigo a julianleioa@gmail.com.